library("knitr") # for knitting RMarkdown
library("kableExtra") # for making nice tables
library("janitor") # for cleaning column names
library("tidybayes") # tidying up results from Bayesian models
library("brms") # Bayesian regression models with Stan
library("patchwork") # for making figure panels
library("gganimate") # for animations
library("GGally") # for pairs plot
library("bayesplot") # for visualization of Bayesian model fits
library("tidyverse") # for wrangling, plotting, etc. Load the poker data set.
df.poker = read_csv("data/poker.csv") %>%
mutate(skill = factor(skill,
levels = 1:2,
labels = c("expert", "average")),
skill = fct_relevel(skill, "average", "expert"),
hand = factor(hand,
levels = 1:3,
labels = c("bad", "neutral", "good")),
limit = factor(limit,
levels = 1:2,
labels = c("fixed", "none")),
participant = 1:n()) %>%
select(participant, everything())Let’s visualize the data first:
df.poker %>%
ggplot(mapping = aes(x = hand,
y = balance,
fill = hand)) +
geom_point(alpha = 0.2,
position = position_jitter(height = 0, width = 0.1)) +
stat_summary(fun.data = "mean_cl_boot",
geom = "linerange",
size = 1) +
stat_summary(fun.y = "mean",
geom = "point",
shape = 21,
size = 4) +
labs(y = "final balance (in Euros)") +
scale_fill_manual(values = c("red", "orange", "green")) +
theme(legend.position = "none")And let’s now fit a simple (frequentist) regression model:
#>
#> Call:
#> lm(formula = balance ~ 1 + hand, data = df.poker)
#>
#> Residuals:
#> Min 1Q Median 3Q Max
#> -12.9264 -2.5902 -0.0115 2.6573 15.2834
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 5.9415 0.4111 14.451 < 2e-16 ***
#> handneutral 4.4051 0.5815 7.576 4.55e-13 ***
#> handgood 7.0849 0.5815 12.185 < 2e-16 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Residual standard error: 4.111 on 297 degrees of freedom
#> Multiple R-squared: 0.3377, Adjusted R-squared: 0.3332
#> F-statistic: 75.7 on 2 and 297 DF, p-value: < 2.2e-16
Now, let’s fit a Bayesian regression model using the brm() function:
fit.brm1 = brm(formula = balance ~ 1 + hand,
data = df.poker,
file = "cache/brm1")
fit.brm1 %>% summary()#> Family: gaussian
#> Links: mu = identity; sigma = identity
#> Formula: balance ~ 1 + hand
#> Data: df.poker (Number of observations: 300)
#> Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
#> total post-warmup samples = 4000
#>
#> Population-Level Effects:
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> Intercept 5.93 0.41 5.12 6.72 1.00 2986 2744
#> handneutral 4.41 0.58 3.30 5.55 1.00 3497 2903
#> handgood 7.10 0.58 5.99 8.29 1.00 3545 2932
#>
#> Family Specific Parameters:
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> sigma 4.12 0.17 3.81 4.46 1.00 3650 2921
#>
#> Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
#> is a crude measure of effective sample size, and Rhat is the potential
#> scale reduction factor on split chains (at convergence, Rhat = 1).
I use the file = argument to save the model’s results so that when I run this code chunk again, the model doesn’t need to be fit again (fitting Bayesian models takes a while …).
Let’s visualize what the posterior for the different parameters looks like. We use the geom_halfeyeh() function from the “tidybayes” package to do so:
fit.brm1 %>%
posterior_samples() %>%
select(-lp__) %>%
gather("variable", "value") %>%
ggplot(data = .,
mapping = aes(y = variable, x = value)) +
geom_halfeyeh()And let’s look at how the samples from the posterior are correlated with each other:
fit.brm1 %>%
posterior_samples() %>%
select(b_Intercept:sigma) %>%
ggpairs(lower = list(continuous = wrap("points", alpha = 0.03)),
upper = list(continuous = wrap("cor", size = 6))) +
theme(panel.grid.major = element_blank(),
text = element_text(size = 12))To compute the MAP (maximum a posteriori probability) estimate and highest density interval, we use the mode_hdi() function that comes with the “tidybayes” package.
fit.brm1 %>%
posterior_samples() %>%
clean_names() %>%
select(starts_with("b_"), sigma) %>%
mode_hdi() %>%
pivot_longer(cols = -c(.width:.interval),
names_to = "index",
values_to = "value") %>%
select(index, value) %>%
mutate(index = ifelse(str_detect(index, fixed(".")), index, str_c(index, ".mode"))) %>%
separate(index, into = c("parameter", "type"), sep = "\\.") %>%
pivot_wider(names_from = type,
values_from = value) %>%
kable(digits = 2) %>%
kable_styling(bootstrap_options = "striped",
full_width = F)| parameter | mode | lower | upper |
|---|---|---|---|
| b_intercept | 6.02 | 5.13 | 6.73 |
| b_handneutral | 4.46 | 3.29 | 5.54 |
| b_handgood | 7.07 | 5.95 | 8.21 |
| sigma | 4.11 | 3.81 | 4.46 |
To check whether the model did a good job capturing the data, we can simulate what future data the Baysian model predicts, now that it has learned from the data we feed into it.
This looks good! The predicted shaped of the data based on samples from the posterior distribution looks very similar to the shape of the actual data.
Let’s make a hypothetical outcome plot that shows what concrete data sets the model would predict:
# generate predictive samples
df.predictive_samples = fit.brm1 %>%
posterior_samples() %>%
clean_names() %>%
select(contains("b_"), sigma) %>%
sample_n(size = 20) %>%
mutate(sample = 1:n()) %>%
group_by(sample) %>%
nest() %>%
mutate(bad = map(data, ~ .$b_intercept + rnorm(100, sd = .$sigma)),
neutral = map(data, ~ .$b_intercept + .$b_handneutral + rnorm(100, sd = .$sigma)),
good = map(data, ~ .$b_intercept + .$b_handgood + rnorm(100, sd = .$sigma))) %>%
unnest(c(bad, neutral, good)) %>%
select(-data)
# plot the results as an animation
p = df.predictive_samples %>%
pivot_longer(cols = -sample,
names_to = "hand",
values_to = "balance") %>%
mutate(hand = factor(hand, levels = c("bad", "neutral", "good"))) %>%
ggplot(mapping = aes(x = hand,
y = balance,
fill = hand)) +
geom_point(alpha = 0.2,
position = position_jitter(height = 0, width = 0.1)) +
stat_summary(fun.data = "mean_cl_boot",
geom = "linerange",
size = 1) +
stat_summary(fun.y = "mean",
geom = "point",
shape = 21,
size = 4) +
labs(y = "final balance (in Euros)") +
scale_fill_manual(values = c("red", "orange", "green")) +
theme(legend.position = "none") +
transition_manual(sample)
animate(p, nframes = 120, width = 800, height = 600, res = 96, type = "cairo")
# anim_save("poker_posterior_predictive.gif")One key advantage of Bayesian over frequentist analysis is that we can test hypothesis in a very flexible manner by directly probing our posterior samples in different ways.
We may ask, for example, what the probability is that the parameter for the difference between a bad hand and a neutral hand (b_handneutral) is greater than 0. Let’s plot the posterior distribution together with the criterion:
fit.brm1 %>%
posterior_samples() %>%
select(b_handneutral) %>%
pivot_longer(cols = everything(),
names_to = "variable",
values_to = "value") %>%
ggplot(data = .,
mapping = aes(y = variable, x = value)) +
geom_halfeyeh() +
geom_vline(xintercept = 0,
color = "red")We see that the posterior is definitely greater than 0.
We can ask many different kinds of questions about the data by doing basic arithmetic on our posterior samples. The hypothesis() function makes this even easier. Here are some examples:
# the probability that the posterior for handneutral is less than 0
hypothesis(fit.brm1,
hypothesis = "handneutral < 0")#> Hypothesis Tests for class b:
#> Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob
#> 1 (handneutral) < 0 4.41 0.58 3.48 5.35 0 0
#> Star
#> 1
#> ---
#> 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
#> '*': For one-sided hypotheses, the posterior probability exceeds 95%;
#> for two-sided hypotheses, the value tested against lies outside the 95%-CI.
#> Posterior probabilities of point hypotheses assume equal prior probabilities.
# the probability that the posterior for handneutral is greater than 4
hypothesis(fit.brm1,
hypothesis = "handneutral > 4")#> Hypothesis Tests for class b:
#> Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
#> 1 (handneutral)-(4) > 0 0.41 0.58 -0.52 1.35 3.05
#> Post.Prob Star
#> 1 0.75
#> ---
#> 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
#> '*': For one-sided hypotheses, the posterior probability exceeds 95%;
#> for two-sided hypotheses, the value tested against lies outside the 95%-CI.
#> Posterior probabilities of point hypotheses assume equal prior probabilities.
# the probability that good hands make twice as much as bad hands
hypothesis(fit.brm1,
hypothesis = "Intercept + handgood > 2 * Intercept")#> Hypothesis Tests for class b:
#> Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
#> 1 (Intercept+handgo... > 0 1.18 0.93 -0.29 2.67 8.5
#> Post.Prob Star
#> 1 0.89
#> ---
#> 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
#> '*': For one-sided hypotheses, the posterior probability exceeds 95%;
#> for two-sided hypotheses, the value tested against lies outside the 95%-CI.
#> Posterior probabilities of point hypotheses assume equal prior probabilities.
# the probability that neutral hands make less than the average of bad and good hands
hypothesis(fit.brm1,
hypothesis = "Intercept + handneutral < (Intercept + Intercept + handgood) / 2")#> Hypothesis Tests for class b:
#> Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
#> 1 (Intercept+handne... < 0 0.86 0.5 0.05 1.68 0.04
#> Post.Prob Star
#> 1 0.04
#> ---
#> 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
#> '*': For one-sided hypotheses, the posterior probability exceeds 95%;
#> for two-sided hypotheses, the value tested against lies outside the 95%-CI.
#> Posterior probabilities of point hypotheses assume equal prior probabilities.
Let’s double check one example, and calculate the result directly based on the posterior samples:
df.hypothesis = fit.brm1 %>%
posterior_samples() %>%
clean_names() %>%
select(starts_with("b_")) %>%
mutate(neutral = b_intercept + b_handneutral,
bad_good_average = (b_intercept + b_intercept + b_handgood)/2,
hypothesis = neutral < bad_good_average)
df.hypothesis %>%
summarize(p = sum(hypothesis)/n())#> p
#> 1 0.03975
Another way of testing hypothesis is via the Bayes factor. Let’s fit the two models we are interested in comparing with each other:
fit.brm2 = brm(formula = balance ~ 1 + hand,
data = df.poker,
save_all_pars = T,
file = "cache/brm2")
fit.brm3 = brm(formula = balance ~ 1 + hand + skill,
data = df.poker,
save_all_pars = T,
file = "cache/brm3")And then compare the models useing the bayes_factor() function:
#> Iteration: 1
#> Iteration: 2
#> Iteration: 3
#> Iteration: 4
#> Iteration: 5
#> Iteration: 1
#> Iteration: 2
#> Iteration: 3
#> Iteration: 4
#> Estimated Bayes factor in favor of bridge1 over bridge2: 3.83284
So far, we have used the defaults that brm() comes with and not bothered about specifiying the priors, etc.
Notice that we didn’t specify any priors in the model. By default, “brms” assigns weakly informative priors to the parameters in the model. We can see what these are by running the following command:
#> prior class coef group resp dpar nlpar bound
#> 1 b
#> 2 b handgood
#> 3 b handneutral
#> 4 student_t(3, 10, 10) Intercept
#> 5 student_t(3, 0, 10) sigma
We can also get information about which priors need to be specified before fitting a model:
#> prior class coef group resp dpar nlpar bound
#> 1 b
#> 2 b handgood
#> 3 b handneutral
#> 4 student_t(3, 10, 10) Intercept
#> 5 student_t(3, 0, 10) sigma
Here is an example for what a more complete model specification could look like:
fit.brm4 = brm(formula = balance ~ 1 + hand,
family = "gaussian",
data = df.poker,
prior = c(prior(normal(0, 10), class = "b", coef = "handgood"),
prior(normal(0, 10), class = "b", coef = "handneutral"),
prior(student_t(3, 3, 10), class = "Intercept"),
prior(student_t(3, 0, 10), class = "sigma")),
inits = list(list(Intercept = 0, sigma = 1, handgood = 5, handneutral = 5),
list(Intercept = -5, sigma = 3, handgood = 2, handneutral = 2),
list(Intercept = 2, sigma = 1, handgood = -1, handneutral = 1),
list(Intercept = 1, sigma = 2, handgood = 2, handneutral = -2)),
iter = 4000,
warmup = 1000,
chains = 4,
file = "cache/brm4",
seed = 1)
fit.brm4 %>% summary()#> Family: gaussian
#> Links: mu = identity; sigma = identity
#> Formula: balance ~ 1 + hand
#> Data: df.poker (Number of observations: 300)
#> Samples: 4 chains, each with iter = 4000; warmup = 1000; thin = 1;
#> total post-warmup samples = 12000
#>
#> Population-Level Effects:
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> Intercept 5.96 0.41 5.15 6.76 1.00 9139 8795
#> handneutral 4.37 0.58 3.23 5.53 1.00 9638 8040
#> handgood 7.05 0.58 5.93 8.19 1.00 9782 8367
#>
#> Family Specific Parameters:
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> sigma 4.13 0.17 3.81 4.49 1.00 12982 9331
#>
#> Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
#> is a crude measure of effective sample size, and Rhat is the potential
#> scale reduction factor on split chains (at convergence, Rhat = 1).
We can also take a look at the Stan code that the brm() function creates:
#> // generated with brms 2.7.0
#> functions {
#> }
#> data {
#> int<lower=1> N; // total number of observations
#> vector[N] Y; // response variable
#> int<lower=1> K; // number of population-level effects
#> matrix[N, K] X; // population-level design matrix
#> int prior_only; // should the likelihood be ignored?
#> }
#> transformed data {
#> int Kc = K - 1;
#> matrix[N, K - 1] Xc; // centered version of X
#> vector[K - 1] means_X; // column means of X before centering
#> for (i in 2:K) {
#> means_X[i - 1] = mean(X[, i]);
#> Xc[, i - 1] = X[, i] - means_X[i - 1];
#> }
#> }
#> parameters {
#> vector[Kc] b; // population-level effects
#> real temp_Intercept; // temporary intercept
#> real<lower=0> sigma; // residual SD
#> }
#> transformed parameters {
#> }
#> model {
#> vector[N] mu = temp_Intercept + Xc * b;
#> // priors including all constants
#> target += normal_lpdf(b[1] | 0, 10);
#> target += normal_lpdf(b[2] | 0, 10);
#> target += student_t_lpdf(temp_Intercept | 3, 3, 10);
#> target += student_t_lpdf(sigma | 3, 0, 10)
#> - 1 * student_t_lccdf(0 | 3, 0, 10);
#> // likelihood including all constants
#> if (!prior_only) {
#> target += normal_lpdf(Y | mu, sigma);
#> }
#> }
#> generated quantities {
#> // actual population-level intercept
#> real b_Intercept = temp_Intercept - dot_product(means_X, b);
#> }
One thing worth noticing: by default, “brms” centers the predictors which makes it easier to assign a default prior over the intercept.
So far, we’ve assumed that the inference has worked out. We can check this by running plot() on our brm object:
Let’s make our own version of a trace plot for one parameter in the model:
fit.brm1 %>%
spread_draws(b_Intercept) %>%
clean_names() %>%
mutate(chain = as.factor(chain)) %>%
ggplot(aes(x = iteration, y = b_intercept, group = chain, color = chain)) +
geom_line()We can also take a look at the auto-correlation plot. Ideally, we want to generate independent samples from the posterior. So we don’t want subsequent samples to be strongly correlated with each other. Let’s take a look:
variables = fit.brm1 %>%
get_variables() %>%
.[1:4]
fit.brm1 %>%
posterior_samples() %>%
mcmc_acf(pars = variables,
lags = 4)Looking good! The autocorrelation should become very small as the lag increases (indicating that we are getting independent samples from the posterior).
Let’s try to fit a model to very little data (just two observations) with extremely uninformative priors:
df.data = tibble(y = c(-1, 1))
fit.brm5 = brm(data = df.data,
family = gaussian,
formula = y ~ 1,
prior = c(prior(uniform(-1e10, 1e10), class = Intercept),
prior(uniform(0, 1e10), class = sigma)),
inits = list(list(Intercept = 0, sigma = 1),
list(Intercept = 0, sigma = 1)),
iter = 4000,
warmup = 1000,
chains = 2,
file = "cache/brm5")Let’s take a look at the posterior distributions of the model parameters:
#> Warning: The model has not converged (some Rhats are > 1.1). Do not analyse the results!
#> We recommend running more iterations and/or setting stronger priors.
#> Warning: There were 1203 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help.
#> See http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
#> Family: gaussian
#> Links: mu = identity; sigma = identity
#> Formula: y ~ 1
#> Data: df.data (Number of observations: 2)
#> Samples: 2 chains, each with iter = 4000; warmup = 1000; thin = 1;
#> total post-warmup samples = 6000
#>
#> Population-Level Effects:
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
#> Intercept 357550121.58 1416057299.71 -2244033111.47 3333594132.43 1.78 3
#> Tail_ESS
#> Intercept 24
#>
#> Family Specific Parameters:
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> sigma 1524412740.64 2392424321.98 21668.93 8317582240.06 1.40 4 41
#>
#> Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
#> is a crude measure of effective sample size, and Rhat is the potential
#> scale reduction factor on split chains (at convergence, Rhat = 1).
Not looking good – The estimates and credible intervals are off the charts. And the effective samples sizes in the chains are very small.
Let’s visualize the trace plots:
fit.brm5 %>%
spread_draws(b_Intercept) %>%
clean_names() %>%
mutate(chain = as.factor(chain)) %>%
ggplot(aes(x = iteration,
y = b_intercept,
group = chain,
color = chain)) +
geom_line()Given that we have so little data in this case, we need to help the model a little bit by providing some slighlty more specific priors.
fit.brm6 = brm(data = df.data,
family = gaussian,
formula = y ~ 1,
prior = c(prior(normal(0, 10), class = Intercept), # more reasonable priors
prior(cauchy(0, 1), class = sigma)),
iter = 4000,
warmup = 1000,
chains = 2,
seed = 1,
file = "cache/brm6")Let’s take a look at the posterior distributions of the model parameters:
#> Warning: There were 3 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help.
#> See http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
#> Family: gaussian
#> Links: mu = identity; sigma = identity
#> Formula: y ~ 1
#> Data: df.data (Number of observations: 2)
#> Samples: 2 chains, each with iter = 4000; warmup = 1000; thin = 1;
#> total post-warmup samples = 6000
#>
#> Population-Level Effects:
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> Intercept -0.13 1.69 -4.10 3.06 1.00 1352 909
#>
#> Family Specific Parameters:
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> sigma 2.04 1.88 0.61 6.95 1.00 1144 1426
#>
#> Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
#> is a crude measure of effective sample size, and Rhat is the potential
#> scale reduction factor on split chains (at convergence, Rhat = 1).
This looks much better. There is still quite a bit of uncertainty in our paremeter estimates, but it has reduced dramatically.
Let’s visualize the trace plots:
fit.brm6 %>%
spread_draws(b_Intercept, sigma) %>%
clean_names() %>%
mutate(chain = as.factor(chain)) %>%
pivot_longer(cols = c(b_intercept, sigma)) %>%
ggplot(aes(x = iteration,
y = value,
group = chain,
color = chain)) +
geom_line() +
facet_wrap(vars(name), ncol = 1)Looking mostly good!
Information about this R session including which version of R was used, and what packages were loaded.
#> R version 3.6.2 (2019-12-12)
#> Platform: x86_64-apple-darwin15.6.0 (64-bit)
#> Running under: macOS Mojave 10.14.6
#>
#> Matrix products: default
#> BLAS: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib
#> LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib
#>
#> locale:
#> [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] forcats_0.4.0 stringr_1.4.0 dplyr_0.8.3 purrr_0.3.3
#> [5] readr_1.3.1 tidyr_1.0.0 tibble_2.1.3 tidyverse_1.3.0
#> [9] bayesplot_1.7.1 GGally_1.4.0 gganimate_1.0.4 ggplot2_3.2.1
#> [13] patchwork_1.0.0 brms_2.10.0 Rcpp_1.0.3 tidybayes_1.1.0
#> [17] janitor_1.2.0 kableExtra_1.1.0 knitr_1.26
#>
#> loaded via a namespace (and not attached):
#> [1] readxl_1.3.1 backports_1.1.5
#> [3] Hmisc_4.3-0 plyr_1.8.5
#> [5] igraph_1.2.4.2 lazyeval_0.2.2
#> [7] splines_3.6.2 svUnit_0.7-12
#> [9] crosstalk_1.0.0 rstantools_2.0.0
#> [11] inline_0.3.15 digest_0.6.25
#> [13] htmltools_0.4.0 rsconnect_0.8.16
#> [15] fansi_0.4.0 checkmate_1.9.4
#> [17] magrittr_1.5 cluster_2.1.0
#> [19] modelr_0.1.5 matrixStats_0.55.0
#> [21] xts_0.11-2 prettyunits_1.1.1
#> [23] jpeg_0.1-8.1 colorspace_1.4-1
#> [25] rvest_0.3.5 haven_2.2.0
#> [27] xfun_0.12 callr_3.4.0
#> [29] crayon_1.3.4 jsonlite_1.6.1
#> [31] lme4_1.1-21 survival_3.1-8
#> [33] zoo_1.8-6 glue_1.3.1
#> [35] gtable_0.3.0 webshot_0.5.2
#> [37] pkgbuild_1.0.6 rstan_2.19.2
#> [39] abind_1.4-5 scales_1.1.0
#> [41] DBI_1.1.0 miniUI_0.1.1.1
#> [43] htmlTable_1.13.3 viridisLite_0.3.0
#> [45] xtable_1.8-4 progress_1.2.2
#> [47] HDInterval_0.2.0 ggstance_0.3.3
#> [49] foreign_0.8-72 Formula_1.2-3
#> [51] stats4_3.6.2 StanHeaders_2.19.0
#> [53] DT_0.11 htmlwidgets_1.5.1
#> [55] httr_1.4.1 threejs_0.3.1
#> [57] arrayhelpers_1.0-20160527 RColorBrewer_1.1-2
#> [59] ellipsis_0.3.0 acepack_1.4.1
#> [61] pkgconfig_2.0.3 reshape_0.8.8
#> [63] loo_2.2.0 farver_2.0.3
#> [65] nnet_7.3-12 dbplyr_1.4.2
#> [67] labeling_0.3 tidyselect_1.0.0
#> [69] rlang_0.4.5 reshape2_1.4.3
#> [71] later_1.0.0 munsell_0.5.0
#> [73] cellranger_1.1.0 tools_3.6.2
#> [75] cli_2.0.0 generics_0.0.2
#> [77] gifski_0.8.6 broom_0.5.3
#> [79] ggridges_0.5.1 evaluate_0.14
#> [81] fastmap_1.0.1 yaml_2.2.1
#> [83] processx_3.4.2 fs_1.3.1
#> [85] nlme_3.1-142 mime_0.9
#> [87] mvnfast_0.2.5 rstanarm_2.19.2
#> [89] xml2_1.2.2 compiler_3.6.2
#> [91] shinythemes_1.1.2 rstudioapi_0.11
#> [93] png_0.1-7 reprex_0.3.0
#> [95] tweenr_1.0.1 stringi_1.4.6
#> [97] highr_0.8 ps_1.3.2
#> [99] Brobdingnag_1.2-6 lattice_0.20-38
#> [101] Matrix_1.2-18 nloptr_1.2.1
#> [103] markdown_1.1 shinyjs_1.0
#> [105] vctrs_0.2.3 pillar_1.4.3
#> [107] lifecycle_0.1.0 bridgesampling_0.7-2
#> [109] data.table_1.12.8 httpuv_1.5.2
#> [111] R6_2.4.1 latticeExtra_0.6-29
#> [113] bookdown_0.16 promises_1.1.0
#> [115] gridExtra_2.3 codetools_0.2-16
#> [117] boot_1.3-23 MASS_7.3-51.4
#> [119] colourpicker_1.0 gtools_3.8.1
#> [121] assertthat_0.2.1 withr_2.1.2
#> [123] shinystan_2.5.0 parallel_3.6.2
#> [125] hms_0.5.3 rpart_4.1-15
#> [127] grid_3.6.2 minqa_1.2.4
#> [129] coda_0.19-3 snakecase_0.11.0
#> [131] rmarkdown_2.0 shiny_1.4.0
#> [133] lubridate_1.7.4 base64enc_0.1-3
#> [135] dygraphs_1.1.1.6